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Abstract. This paper presents a detailed asymptotic and numerical investigation of the phase 
diagram for global minimizers to a Cahn-Hilliard functional with long-range interactions in two 
space dimensions. We introduce a small parameter measuring perturbation from the minimal order- 
disorder transition, and derive asymptotic estimates for stability regions as the parameter tends 
to zero. Based upon the gradient flow, we introduce a hybrid numerical method to navigate 

through the complex energy landscape and access the ground state of the functional. We use this 
method to numerically compute the phase diagram. Our asymptotic predictions show surprisingly 
good agreement with our numerical results. 
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1. Introduction. In this article, we asymptotically and numerically address 
the phase diagram with respect to the parameters 7 and m for the following mass- 
constrained variational problem: For 7 > and m G (—1, 1), minimize 



(^|V^i|2 + iL^^)dx -I- / / G{x,y)iuix)~m)iuiy)~m) dxdy, (1.1) 



0^7^' 4 / JnJn 

over all u with j^^-^ udx = m. Here G denotes the Green's function of —A on a cubic 
domain 57 :— [0, L]" C M" with periodic boundary conditions. We refer to functional 
(1.1) as a Cahn-Hilliard functional (c/. [1]) with long-range interactiontj^ It may 



simply be viewed as a mathematical paradigm for energy-driven pattern formation 
induced by competing short and long-range interactions: Minimization of the first 
two terms (short-range) leads to domains of pure phases of m = ±1 with minimal 
transition regions, whereas the third (long-range) terms induces oscillations between 
the phases according to the set volume fraction m. On a sufficiently large domain 
f2, the competition of the two leads to pattern formation on an intrinsic scale which 
depends entirely on 7. Throughout this article, we always take the physical domain 
H. of siz(Qmuch larger than this intrinsic scale (c/. Figure 1.1). 



A tool for our analysis is the H ^ gradient flow (c/. ^) for u := u — m which 
takes the form 

CJU 1 

— = A^ u A (li^ + imu^ - (1 - im^)u) - u. (1.2) 

ot ^■^ 

It is important to note here that we compute the gradient with respect to H~^, a 



nonlocal metric. Hence the presence of the nonlocal term in the functional ( 1.1 ) simply 
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Fig. 1.1. Typical long-time solutions to fLBy with 7 = 10. Left: Lamellae, m = 0. This is 
a global minimizing state. Center: Hexagonally packed spots, m = .3. This is a global minimizing 
state. Right: A mixed state, m = .4. This is a typical metastable solution. All figures were computed 
on a domain of size 4it X An. Here and in the following figures u = 1 is represented in black and 
u = —I in white. 



gives rise to a local perturbatiorj^ of the standard Cahn-Hilliard equation. However, 
it significantly changes its behavior, making it in some ways a hybrid of the Swift- 
Hohenberg equation |Mj and the standard Cahn-HiUiard equation. 



Whilst the functional (1.1) is mathematically interesting on its own, we were 
drawn to it because of its direct connection with self-assembly of diblock copoly- 
mers: the functional is a rescaled version of a functional introduced by Ohta and 
Kawasaki jJH [23] . Melts of diblock copolymer display a rich class of self-assembly 
nano-structures from lamellae, spheres, cylindrical tubes, to double-gyroids and other 
more complex structures (see for example, [HHl]). Moreover, the usefulness of block 
copolymer melts is exactly this remarkable ability for self-assembly into particular 
geometries. For example, this property can be exploited to create materials with de- 
signer mechanical, optical, and magnetic properties [2]. Therefore from a theoretical 
point of view, one of the main challenges is to predict the phase geometry /morphology 
for a given set of material parameters, that is, the creation of a phase diagram. 

The state of the art for predicting the phase diagram is via the self-consistent mean 



field theory (SCFT) PUIIII]. While the simple functional ( 1.1 ) can be connected with 
the SCFT via approximations (c/. [TU]) with increasing validity close to the order- 
disorder transition, it is generally regarded as the basis of a qualitative theory, and 
one might question its usefulness with regard to predicting self-assembly structures 
for given material parameters. Preliminary numerical experiments [351 [HI indi- 
cate that all the phases (including double gyroids and perforated lamellae), some of 
which had been predicted using the SCFT [5D] and all of which have been observed 
for polystyrene-isoprene '15|, can be simulated as minimizers of ( |1.1[ ) starting from 
random initial conditions. This begs the question as to the extent to which a phase 



diagram via ( 1.1 ) can be compared with those of experimental observations and SCFT 
calculations, at least close to the order-disorder transition. 

A thorough numerical phase diagram for ( |1.1[ ) with n — i (3D) is by no means an 



easy task. In addition to numerical complications associated with the stiff PDE ( 1.2 ) 



and the necessary small time steps and large spatial grids, the energy landscape of 



( 1.1 ) is highly non-convex, with multitudes of local minimizers and metastable states 



about which the gradient flow dynamics are very slow. Many of the the simulations we 



presented in l^j, were not simply the final steady states for simulations of (1.2) with 
random initial conditions. Rather, one tended to get stuck in metastable states and 
some procedure for exiting these metastable states in order to flow to lower energy 



•^Note that the gradient of the nonlocal term with respect to L^ would be (—A) ^(«), and 
working in has the effect of introducing an additional (—A). 
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states was crucial. Such procedures are often loosely dubbed simulated annealing. 
The present article may be viewed as a test case for the 3D phase diagram. Here, 
via a combination of asymptotic analysis and numerical experiments, we consider the 
phase diagram in space dimension n = 2 with respect to both global and stable local 
minimizers. The 2D situation is greatly simplified as the range of possible minimiz- 
ers is both drastically reduced and, in fact, well accepted to be only certain basic 
structures: 



• disordered, i.e. minimizers of (1.1 ) are simply the uniform state u = m; 

• lamellar, i.e. minimizers of (1.1) have a one dimensional structure; 

• spots, i.e. minimizers of (1.1) are a periodic array of (approximate) circles 
arranged on either a hexagonal or rectangular lattice (in what follows, all 
spot solutions are hexagonally packed unless otherwise indicated). 

Examples of the latter two cases and a metastable mixed state are presented in Figure 

We adopt periodic boundary conditions throughout this article and view the do- 
main as a flat torus. Our numerical and asymptotic results demonstrate explicit 
regions in the 7 vs m plane for stability of the lamellar, circular and disordered states 
and explicit regions in the 7 vs m plane wherein the lamellar, hexagonally packed 
circular, and disordered states are global minimizers. At the surface, this may not to 
appear to be a particularly novel contribution. Indeed, there are many PDE models 
for pattern formation, some of which, like this one, are variational, and transitions be- 
tween disorder, spots and lamellae as volume fraction changes and asymptotic-based 
descriptions are not surprisingly ubiquitous (c./. [251 IS])- For instance, there is 
considerable interest in localized solutions to the Swift-Hohenberg equation and other 
similar equations in dimensions 1, 2 and 3 - see for example 34, 3, 31 , 18. i25lfTBl[?flfT^ 
and the references therein. In periodic solutions to the Swift-Hohenberg equation 
are analyzed in one dimension with both the detuning parameter a and the imposed 
fundamental period length as parameters. For our problem, the period length is de- 
termined by the minimization procedure, and not solely by a parameter. More impor- 
tantly, our focus here lies on the two-parameter, highly non-convex energy landscape 



of the general functional (1.1 1, not simply on steady state solutions to the partial 



differential equation (1.2 1. Specifically, the novelty of our approach resides in the: 



(i) the demonstration of a hybrid numerical method to address the complex 



metastability issues and provide access to the ground state of (1.1 ); 
(ii) the (seemingly surprising) strong agreement (c.f. Section [5| of our numerical 
results with the results of standard asymptotic analysis suggesting that for the 
purposes of describing the basic geometric morphology of the ground state of 



(1.1), the linear behavior of (1.2) gives very accurate predictions in a (finite) 

neighborhood of the order-disorder transition 
To this end, we are not aware of any systematic study for similar models with two 
parameters in two or more dimensions. Finally, we note that both (i) and (ii) above 
are particularly important with regard to a 3D study wherein the set of possible 
candidates for minimizers (global and local) is far more complex and in fact unknown. 

2. Asymptotic calculations. In this Section we identify local and global min- 



imizers of the energy functional ( 1.1 ) by determining the stationary solutions to ( 1.2 1 
and evaluating their stability and energy. Asymptotic methods are used to describe 
the solutions in the limit [7 — 2| <|; 1 and numerical experiments show they are valid 
over a much larger parameter region. We begin by identifying the order-disorder tran- 
sition (ODT) curve which separates regions of phase space where u = m is linearly 
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Fig. 2.1. Asymptotic coordinates. We identify points in the plane via (firn* {■y) , -y) where 
m*(7) <C 1 as 7 4, 2. 



stable from those regions where it is not. 

Consider the hnearisation of (1.2) about u = Q 

vt^Cv = — \rA'^v-(l-3m'^)Av-v (2.1) 
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where v is periodic on Lx/2] x [—Ly/2, Ly/2]. Taking the Fourier transform 

shows that the eigenvalues of C are given by: 



A = -^\ik,, ky)\^ + (1 - 3m^)\{k^,ky)\^ 



where (k^, ky) represents the two-dimensional wavevector. The condition max(Re(A)) 
determines the ODT curve as 



7-2 

1 — 3m^ V 87 



7 = or 



At fixed 7, u = is linearly stable for m > m*{-^) and unstable for m < m*(7). On 
the ODT curve there is a finite dimensional center manifold on which we will construct 
solutions and determine their stability. For /3 > 0, we set m = /3m* (7) and identify 
points in the plane via (771,7) where m*(7) <C 1 as 7 J, 2 (see Figure 2.1). Another 
of way of saying this is that we consider the asymptotics as (m,7) approaches the 

le form 7 = — ^ 

expansion for u for small to*: 



point (0, 2) along curves of the form 7 = — ^—2-. We introduce a regular asymptotic 

1-3^ 



u 



{x,y) m*ui{x,y) + {m*Yu2{x,y) + . . . , (2.2) 



where ui ^ 0, since we are interested in the nonzero solutions of (1.2). Consideration 
of the ©(to*) equation along with periodic boundary conditions leads to the following 
leading order representation of u for to* ^ 1, 

u{x,y,t) ^ m*{a{t)(pi{x,y) + h{t)(j)2{x,y) + c{t)(l)3{x,y)) , (2.3) 

where 

<j>i — a/2 cos (^V2x^ , 02 = V^cos ( — -^a; 





Fig. 2.2. Asymptotic results. Stability diagram of stationary solutions (left). Solid lines rep- 
resent global stability boundaries and dashed are linear .stability boundaries. Lamellae are globally 
stable between m = and the first solid line and linearly stable until the dashed line. Spots are 
globally stable between the solid lines and linearly stable between the dash-dotted lines. The dotted 
line denotes the linear stability transition of the constant state. Energy \2.i^ for steady solutions of 
(right). 



This expansion captures all the possible symmetric periodic configurations described 
above. We proceed to determine the amplitude dynamics on the center manifold of 
( plj ) X" = span{0i,(/)2,(/)3} by projecting the fuU PDE onto Thus, we 

consider 

(^ut + \ A^u - A (u^ + 3l3mu'^) + (1 - S/^^m^) Au + u, (j)^ = 0, 

for i = 1,2,3 on the domain i7, appropriately defined by = ^ and Ly = -^L^. 
Computing the inner products in and expanding in powers of m* with the aid of 
Maple, we find the following amplitude ODE system 

a = 6(1 - I3'^)a - 6V2pbc - 6(6^ + c^)a - 3a^ 

& = 6(1 - (3^)b - 6V2;3ac - 6(0^ + c^)b ~ 36^ (2.4) 



c = 6(1 - P^)c - 6V2l3ab - 6(0^ + b^)c - 3c^ 



where time has been rescaled with (m*)^. 

The original equation is a gradient system, and thus, the reduced system is one 
as well. The corresponding Lyapunov function is given by 

V{a,b,c) = - I3^)(^a^ + b^ + c^) + 6\/2/3a6c 

+ 3{a^b^ + b^c^ + a^c^) + ha^ + b^ + c^). (2.5) 



This function is in fact the projection of the energy functional (1.1) onto the center 
manifold X'^. We consider the structure and stability of the stationary solutions of 
the amplitude ODE system by analyzing its Lyapunov function. The stationary states 
satisfy the following system 

Va{a,b,c) = Q, V6(a,6,c) = 0, and Vc(a, &, c) = 0, (2.6) 

and they are linearly stable when all eigenvalues of the Hessian matrix, H{V{a, b, c)), 
are positive. We identify five fixed points of the system (2.6 1: 
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1. a = 0, 6 = 0, c = — 7- disorder; 

2. a = ±^2(1-/32), 6 = 0, c = lamellae; 

3. a ^ b = ±y^^^^-^^-^, c = — > Triangular spots; 

4. a = a, 6 = c = ±a, where a ~ 
circular; 



P ± -^/S — 4/32 — )■ hexagonally packed 



5. a = 6 = ±a, c = — c, or a = — 6 = ±a, c = c, where a = y^ 2(i ^sff ) 
c = 2^/2(3 — |a| = 6 7^ c case. 
While other similar systems can generate other patterns (c/. [H]), none of these 
appear in this system. Evaluating the eigenvalues of H{V{a,h,c)) at the five fixed 
points, we determine the following three regions of linear stability: 

1. < /3 < lamellae; 

2. -^j^ < /3 < hexagonally packed circular, 

3. $ > I: disorder. 

Evaluation of the Lyapunov function at the three linearly stable steady states is then 
used to determine the following three regions of global stability: 

1. < /? < jgV551 - 174\/6: lamellae; 

2. 27)\/551 — 174\/6 < 13 < 3 ^/^^: hexagonally packed circular, 

3. /? > 3-y/^: disorder. 

Note that there are regions where more than one solution type is linearly stable. 
Figure [2?2| shows the asymptotic linear stability and global stability diagrams. 

3. Numerical Simulations. In this section, we present the results of our nu- 
merical experiments and compare them with the asymptotic calculations of the previ- 
ous section in regimes where they may present some validity. Our numerical method 
will be described in some detail in Section |4] It is a hybrid method which not only 
integrates the PDE (1.2) but also involves an interplay with the energy through 



certain methods of simulated annealing. For most values of the parameters (771,7), 



believe that this method results assesses the ground state of (1.1), at least from the 
point of the it's inherent structure (geometry and symmetry). Of course we have no 
proof of this statement, and moreover, rigorous results pertaining to the ground state 
are difficult to obtain even in 2D (see [HllHlElllSDj for partial results). 

We sampled the parameter plane by taking a sequence of randomly chosen points 
(mi,7i) € [0,1] X [2,25] and implemented our method for each such (771^,7.^) with 
u{xi,yj,t = 0) G (—1,1) randomly chosen from a uniform distribution. Once a 
sequence of more than 500 runs was complete, an edge detection algorithm was used 
to place additional points near the interfaces between regions where u = 0, hex spots 
and lamellae are stable. For all runs we computed the energy ( |1.1[ ) at each time step, 
making sure that the final presented state had the least energy over the course of the 
entire run. 

3.1. Results of numerical simulations. Our numerical experiments confirm 
the asymptotic description of the stability curves and the solution structure in the 



limit 7^2, to* (7) | 0. In Figure 3.1 we present the numerically computed bifurcation 
diagram with the asymptotic stability curves overlaid. Here x , o and o indicate stripes, 
hex packed spots and the disordered state respectively. The solid lines mark the 
global stability curves and the dashed lines outline the linear stability regimes. Note 
that almost no runs converged to a state outside of its asymptotic region of global 
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Fig. 3.1. (These figures need to be viewed in colour) Numerically computed phase diagram. 
(Bottom) Complete diagram. (Top) Detail for ^ close to 2. Blue crosses: Lamellae. Red circles: 
Hex packed spots. Black diamonds: disorder. The red dashed-dotted lines mark the linear stability 
boundary of spots, the blue dashed-dotted line marks the linear stability boundary of lamellae, the 
black dashed-dotted line marks the linear stability boundary of the disordered sate, and the solid 
black lines mark the global stability regions of lamellae and spots respectively. 
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Fig. 3.2. Comparison of energies along solution branches for 7 = 2.5, 3, 3.5, 5, 10 and 20. Note 
that the continuation routine failed for some values of 'y on the stripe branch for m sufficiently large 
and on the spot branch for m sufficiently small but these always occurred well beyond the value of m 
at which the branches exchange global stability. In each energy diagram, the lower value of 7 has the 
higher energy. The vertical bar is the asymptotic value m* (7) where the branches exchange stability. 
On the right, there are three plots for each run from random initial data: m < m* , m ~ m* , m > m* . 
Here we can see that there is no convergence sufficiently close to m* but the expected patterns appear 
away from it. Notice that the there is greater separation in the colors and hence the solution extremes 
as 7 increases. 
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stability and no runs converged to a state outside of its region of linear stability. 



Figure 3.1 (bottom) displays a detailed examination of the bifurcation diagram for 
2 < 7 < 5. This is within the expected range of asymptotic validity and shows very 
good agreement. 

Figure |3.1| (left) has 2 < 7 < 25 and the asymptotic estimates are now only 
in qualitative agreement as 7 increases. Limited additional runs up to 7 = 100 
were also performed and still show qualitative agreement with asymptotics in that no 
new phases were seen and the ordering of stripes, spots and the homogeneous state 
occur for increasing m. Unfortunately, the stability region of the PDF time-stepper 
quickly decreases with increasing 7 and the time-scale of PDF evolution slows down 
making the runs with very large 7 increase in cost faster than 7^ and hence very 
computationally expensive. For increasing 7 we find that the stripe/spot transition is 
much better approximated by the asymptotics than the spot/homogenous transition. 
That is, for increasing 7 there is an ever wider region where spots are globally stable 
but the homogeneous state is linearly stable. 

Continuation in m of spots and stripes was also performed for 7 — 2.001, 2.01, 2.1, 
2.25, 2.5, 3, 3.5, 5, 10 and 20. Comparing the energies of the different phases allows 
us to identify the global minimizer directly. Figure [3?2| shows that the agreement with 
the direct PDF computations is very good. However, computationally continuation 
is much slower and is not completely reliable. Continuation keeps the geometry fixed 
as m varies not necessarily guaranteeing a global minimizer over the whole of the 
branch. Typically, we computed three branches of solutions and kept the ones with 
lowest energy in the region of the transition. 

These results clearly demonstrate that the asymptotic analysis presented in Sec- 
tion 2 accurately describes both the phase diagram and solution structure of minimiz- 
ers beyond the point (to, 7) — (0,2) and up to perhaps as far as 7 = 5. For 7 > 10, 
the agreement is much weaker between the asymptotics and numerics. 

4. Numerical methods. The details of many aspects of our numerical imple- 
mentation are straightforward and well understood. We use a pseudo-spectral method 
in space, because of the large scale periodicity, and two different timestepping schemes 
for evolution. For early times, we use exponential time-differencing (FTD) as the dy- 
namics are quick and we need small timesteps to resolve them. FTD provides a cheap 
per step highly accurate method with well known stability and accuracy properties 
for stiff diagonalizable PDFs. For later times we switch to an iterative linearly im- 
plicit gradient stable algorithm [37] . This method is more expensive per step and only 
first order accurate but it allows arbitrarily large time steps and guarantees that the 
energy will not increase. 

For all computations we used N — 256 spatial modes in both dimensions and 



= 1+^3/2 ■ The initial choice of domain size was somewhat arbitrary (c/. (4.1 1). 
As we have mentioned the domain size used is sufficiently large enough with respect 
to the intrinsic period scale which is determined by 7. However, even with such a 
choice, the exact size of the domain can influence both the minimizing patterns and 
certainly the gradient flow route towards them. Thus we have accounted for this in our 
algorithm via a variation of domain size. This is discussed in Section |4.4[ Typically 
we computed for < t < 100 = tp but took tp a.s large as 2500 in some cases. All 
computations were done on a laptop and implemented in Matlab. Unless otherwise 
indicated, the norm of the residual (||ut||2) is less than 10^* in all figures. 
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Fig. 4.1. Metastable state with 7 = 20,m = 0.8. Left: t = 100. Right: t = 10000. The 
difference in energies between these profiles is less than le — 7 but clearly the left state is not an 
energy minimizer as the packing of spots is irregular. The dynamics is driven by weak interaction 
between the spots. 




Fig. 4.2. Comparison of two runs with the same initial conditions. The top simply integrates 
the PDE while the bottom implements the spectral weighting algorithm described in the text. Both 
runs are identical for <t < 40. At t = 40 the spectral weighting is turned on and acts very quickly 
to allow lamellae to develop. Note that at this scale the difference in the energies is not noticeable 
(see Figure \4.3\below). In both runs 7 = 10 and m = 0.1. 




20 40 60 80 100 120 



time 

Fig. 4.3. Detail of energy over time for Figure \4-.^ The dashed line has no spectral weighting. 
The solid lines all have spectral weighting for 40 < t < 80. The damping parameter p takes the 
values .05, .1, .15, .2, .25 and .3. Over this range this is no discernible difference in the final states 
with weighting. For all runs 7 = 10 and m = 0.1. 
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4.1. Metastability. Figure [42| shows snapshots of a typical run as well as the 
energy decay over time. Notice that the "final" profile is a labyrinthine-like pattern 
rather than pure lamellae or spots but that the energy decay is, relatively, small at 



that time. This phenomenon is pervasive for this problem, as the functional (l.ll 
has many local minimizers and metastable states near which the dynamics are very 
slow. Numerically, one cannot distinguish a metastable state from a stable one, since 
they are both identified as solutions for which the relative change in m or i? between 
timesteps is smaller than some tolerance level. Furthermore, irregular long lasting 
states are common in diblock copolymer experiments I14j , so, without additional 
analysis, it is unclear whether they are steady states, rather than just persistent 
intermediate profiles. 

Analytical considerations can provide insight into the classification of observed 
profiles into metastable and stable steady states. The analytical results on the energy 



functional (1.1) are suggestive that its global minimizers are nearly periodic profiles 
with level surfaces with nearly constant mean curvature T] [301 dU [51 15] . The 
irregular states that we observe do not match this description, and thus we term them 
as metastable states. Long-lived transient behavior due to the metastability of the 
Cahn-Hilliard equation is well-known (see for example, |33j ) . and the metastability is 
also observed in SCFT numerical calculations [13]. Given that our PDF differs from 
the Cahn-Hilliard equation only in the extra nonlocal term, and can be connected to 
the SCFT theory (albeit, via further approximations), there is little reason to believe 
that it can escape the issue of metastability in general. Hence, we need to modify the 
gradient descent to better approximate the true global minimizers. 

4.2. Selecting and damping modes. Techniques for dealing with metastabil- 
ity and highly non-convex energy landscapes often belong to the broad class of statis- 
tical methods, called simulated annealing. They were created to navigate through a 
complex energy landscape in search of a global minimizer. A very simple form of sim- 
ulated annealing can be achieved by adding unbiased noise to the evolved metastable 
state. This may force the solution out of the local minimizer that it is stuck in and 
make it continue its evolution through the energy landscape. Unfortunately, this ap- 
proach does not provide a guaranteed way of addressing metastability as: too much 
noise leads to the divergence of the solution and even when the solution remains 
bounded, there is no way of ensuring that it will not revisit the local minimizers that 
it was stuck in before. Also, the added noise is very quickly damped out due to the 
fourth-order derivative term in this problem. 

A different approach to the removal of defects is provided by the technique of 
spectral filtering. The essence of this method lies in the removal of insignificant spec- 
tral components from an evolved state. That is, we evolve the solution from random 
initial conditions until a structure is formed, compute its Fourier coefficients, and 
keep only the modes which correspond to the coefficients above a certain threshold. 
The evolution is then continued and the process of spectral filtering repeated. This 
approach was suggested in [Hj, and, combined with adding noise, applied to our PDF 
in a modified form. 



Functional ( 1.1 1 can be thought of as a length selection mechanism and, in fact, if 



we plot the energy concentration in Fourier space over time we see that this happens 



quite quickly. Figure 4.4 (left) shows v{k*,t) such that maxk v{k,t) — v{k*,t). In 
Figure 4.4 (center) we see that the energy is concentrated in a narrow band about k* 
in Fourier space. Motivated by this observation we evolve the original PDF until past 
the point that the dominant length-scale has emerged. Denoting this length scale k* 
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Fig. 4.4. (left) k* over time. Here we see that the evolution quickly converges to a dominate 
length scale and that the spectral weighting does not change it. (middle) before the application of the 
spectral weight most energy is concentrated near the dominant mode k* . (right) The effect of the 
spectral weight is to allow the energy to concentrate even more so on k* and its integer multiples. 



we damp the Fourier coefficients as follows: v[k) — > w{k] k*)v{k) where 

w{k;k*) = (1 -p) + 

p (exp (-5(1 - \k\/k*y) + exp (-5(2 - \k\/k* f) + exp (-5(3 - \k\/k* f)) 

This keeps information at all wavelengths but focusses the dynamics at the key length- 
scale and its higher harmonics. Experimentation with the parameter p indicates that 
there is little difference in the outcome with 0.05 < p < 0.3 as indicated in Figure 
|4.3[ With p too small there is no effect and the wrong pattern may emerge if p is too 
large or the standard gradient flow is not run long enough. 

This approach allows us to first identify the energy minimizing length scale, the 
dynamics to focus on these wave lengths, ensure the local stability of profile which 



emerged and then finally to smooth out the effects of added noise. Figure 4.2 shows 
two runs from the same initial conditions. The top has no spectral damping and ends 
up stuck in a metastable state whereas the bottom leads to the global minimizer. The 
energy of both runs is shown in the middle Figure |4.2| where the difference between 
the two runs is almost indistinguishable. A detailed view of the energy is presented 



in Figure 4.3 Lastly, Figure 4.4 (right) shows the coefficient distributions of the final 
time profiles in Figure [4^ with + for global minimizing state and • for the metastable 
one. Notice that the spectral weighting did not shift k* but rather allowed a simpler 
pattern to emerge. 



Remark We have also included four movie files that clearly demonstrate the use of 
this approach. We compare two cases (to, 7) = (.1, 10) and (to, 7) = (.25, 10) showing 
the effect of taking p = or p = .1. For to = .1, the unmodified run leads to a 
metastable mixed state including both spots and stripes whereas modified run leads to 
pure stripes with lower energy. Form — .25 the unmodified run leads to a non-optimal 
collection of spots of differing sizes distributed seemingly randomly. Here, the modified 
run leads to uniform spots packed on a hexagonal grid. All four runs were started from 
the same random initial data. 

4.3. Choice of time-stepping method. When deciding on the most appro- 
priate time discretization scheme for numerical simulations of the PDF, we considered 
exponential time-differencing methods of different order and construction as well as 
two gradient stable schemes based on the work of Eyre [T3] . We settled on the ET- 
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DRK4 scheme [17 as it has the highest order of accuracy and its stabihty region is 
larger than that of the ETD2 scheme. However, the time step is stabihty hmited so to 
compute for long times when the dynamics are slow we need something else. For large 
timesteps we found that a fully implicit gradient stable method |131 137j was the most 
robust in general and the fastest when we could take large timesteps. Combining the 
two methods led to an algorithm at least 3 times faster than either one alone. 

The gradient stable scheme is stable for all timesteps and guarantees decay of the 
energy but requires a nonlinear solve at each timestep. The authors in [37] developed 
an iterative scheme to do this but the number of iterations required greatly depends on 
the magnitude of the solution change over each timestep. This means it is not feasible 
to use it to take very large steps at the beginning. To further speed the computation 
we developed an adaptive timestep routine based on the number of required iterations. 

We also considered numerical implementation of the volume constrained gra- 



dient flow of the energy functional ( l.ll. This consideration was motivated by the fact 



that this evolution system is less stiff than the PDE (1.2 1, as it contains a laplacian, 
rather than a bilaplacian term. However, we found that it did not result in a more 
efficient computation of the steady states as it required far more computational steps. 
This could be because the lower order derivatives do not damp out high frequency 
modes as quickly. 

Continuation in m was done with a modified Newton's method allowing for lin- 
early implicit iterations without computing the 256^ x 256^ Jacobian matrix. 

4.4. Domain size selection. We do not know the natural period length of 
arbitrary solutions and we do not want to enforce a solution type due to our choice of 
box size thus we adaptively choose the domain size. In the limit 7 J, 2 the asymptotic 
form of the solution has |fc| = 1/2 so we take Lq — 12it as a preliminary box size to fit 
a reasonable number of periods in the box. As 7 increases, the intrinsic period will 
decrease. For large 7, this intrinsic length scales like (1/7)^/'^ (c/. [71 [301 11] ) and 
hence as a very rough measure, we found it useful to set the initial domain size as 

2x1/3 



= Lo (^-j . (4.1) 

However, we want to ensure that we have reached a true global minimizer, and that 
boundary effects resulting from a "bad" choice of domain size have not prevented us 
from reaching this goal. Thus at the end of each run, we consider the field u which has 
been numerically computed on [0, L]^, rescale it to a unit domain [0, 1]^, and examine 
the rescaled form of the energy per unit area: 



^[o,i]2 7^ ■ J[Q.,iY 

/ G{x,y) {u{x) — m) {u{y) — m) dx dy 

i[0,l]^ "'[0,1]^ 

Ii{u) + hiu) + L^h{u,m). 



We then minimize E with respect to L for fixed £t, 7 and m. With this choice of 
optimal L, we then integrate for a small time further to ensure that the energy truly 
is lower. This also gives us an additional mechanism to find meta-stable states as 
they will have either a higher energy per unit volume or drastically different optimal 
box size than their near neighbors in the (m, 7) plane. 
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4.5. Automatic pattern identification. To aid in the automatic classification 
of evolved states we examine the distribution of peaks with similar wave length to that 
of the dominant mode. In particular, we define and analyze the following function 



m = / exp 



\{kx,ky) - (sin0,cos0)fc* 



-v{k) dk. 



If g{9) (where 9 = tan ^ |^ is the polar angle) is found to have two dominant peaks 
for < < 27r, the state is classified as a stripe. If it has six then it is classified as a 
hexagonally packed circular state. By symmetry, these are the n = 1 (lamellae) and 3 
(hex packed spots) cases from Section 2 respectively. Labyrinthine structures do not 
display either of these patterns. This method is not perfect, working only in about 
95% of cases. 



4.6. Summary of algorithm. In summary, we integrate ( 1.2 ) using the follow- 
ing algorithm: 

1. Set L = f I j 127r. 

2. Choose random initial data in [— 1 — m, 1 — m] with mean m. 

3. Integrate Ol) using ETDRK4 until t ^ h. 

4. Integrate (1.2) using ETDRK4 with spectral weighting for ti < t < t2- 

5. Integrate (1.2) using ETDRK4 with added white noise for t2 < t < t^. 

6. Integrate (1.2) using the nonlinear gradient stable scheme for < t < t^. 



Domain size variation/selection of Section 4.4 



9. Classify the final state. 
This procedure may appear overly complicated but recall that we are attempting 
to minimize a non-convex, non-local energy in a manner which does not use an antici- 
pated solution structure (as this is not possible in 3D), and in a way which is efficient, 
reliable and automatic. To do this we have combined two different time-stepping 
strategies with two different annealing mechanisms. 

5. Unusual Agreement with Asymptotic Analysis. The calculations of the 
global stability curves of Section [2] give surprisingly good agreement with the results 
of our spectrally weighted numerical algorithm for a larger range of 7 than might be 
initially expected. This agreement extends far beyond the immediate neighborhood of 
the order-disorder curve. Naturally, this should raise a certain amount of skepticism 
and concern that the weighted numerical method biases the solutions towards the 
linear dynamics. Let us address these concerns. First, we note that the spectral 
weighting is only employed for part of run. For those parts, it is certainly possible 
that a certain bias is introduced towards the linear dynamics. However, recall that 



our fundamental goal here is the phase diagram of ( 1.1 ): that is in parameter space 7 



vs m, characterizing the geometric morphology of the global minimizer. Even without 
the spectral weighting, all runs gave final states in which the basic pattern dichotomy 
(spots vs stripes) was recognizable. The spectral weighting allowed for final states 
with lower energy and which were identical to the desired morphology (eg. straight 
stripes vs corrugated). At certain stages of the evolution, it does introduce a bias 
towards lowering the energy and if this bias favors the "linear regime" of (1.2 1, so 
be it. In fact, all of these conclusions lead to the following observation: For the 



purposes of describing the basic geometric morphology of the ground state of (1.1) 



the linear behavior of (1.2) gives very accurate predictions in a (finite) neighborhood 
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7 


m 


mecLn^(max(ii) — min(ti)) 
2 


2.001 


0.0050 


0.013 


2.01 


0.0157 


0.043 


2.1 


.0485 


0.109 


2.25 


0.0741 


0.163 


2.5 


0.2582 


0.240 


3 


0.3333 


0.317 


3.5 


0.3780 


0.413 


5 


0.4472 


0.626 


10 


0.5164 


0.956 


20 


0.5477 


1.021 



Table 5.1 

Variations about the mean as m* varies along the energy minimizing branches (see Figure \3. S\ 



of the order-disorder transition. This is particularly true for initial data which is 
a small perturbation of the homogeneous state. That is, the dynamics are initially 
driven by the linear terms and then there is a large energy barrier to overcome to 
transition from one phase to the other. In fact, starting with small initial data gives 
convergence to states as predicted by the asymptotics for 7 as large as 25. Starting 
with initial conditions with large random deviations generates the diagram presented 
here as it is with these less biased initial conditions that we our able to access the lowest 
energy configurations. The domination of the linear terms in constructing stationary 
solutions to higher-order PDE is well known in the folklore but little studied beyond 
formal observations [S] 

However, there is another more fundamental reason why this problem is unusual: 
The energy has a term which keeps solutions in the linear regime. Solutions to the 
standard Cahn-Hilliard equation immediately tend to values of ±1 -|- C(£) where e is 



the interfacial thickness. In (1.1 ), the long-range interactions prefer to keep solutions 
close to oscillations about the mean. But this is precisely the asymptotic Ansatz, and 
thus it is no surprise that it does well to capture the solutions when they are still 
in this regime. Table [5T] presents the magnitude of fluctuation about the mean as 7 
varies and shows there is a large range of 7 where the solution is in the asymptotic 
regime. 

There is also a region of the phase diagram where the agreement is quite poor: 
for large 7 and large m there is a region where there are very well separated hex 
packed spots and a very small region where there are Cartesian packed spots as well. 



Examples of both of these with 7 = 20 are presented in Figure 5.1 Neither of these 
behaviors appears for small m* and the Cartesian packed spots are predicted to never 
be local minimizers. However, SCFT does predict a region of "close-packed" spheres 
which would be the 3D analogue of this behavior. 

6. The Outlook for 3D. The 3D problem is considerably more complicated 
because of the multitude of both local and global minimizers. Preliminary experi- 
ments suggest that our numerical method is reasonably successful throughout much 
of the phase plane. Phase boundaries on the other hand are much more difficult to 
capture and it is thus very encouraging that, at least in 2D, asymptotic analysis pro- 
vides such a good indication of where these are close to the oder-disorder transition. 
Naturally, calculations for double gyroid and perforated lamellar symmetries are more 
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Fig. 5.1. Unusual solutions not predicted by linear analysis. Both these cases have energy lower 
than the mixed state but are not predicted by the linear analysis. 



complicated but are still tractable. As we saw in [5] (though the simulations there 
were performed with a simpler, more ad hoc numerical method), metastability issues 
in 3D are even more complex and additional methods of simulated annealing may 
invariably be needed for a full phase diagram close to the order-disorder transition. 
Moreover our blind test for characterization of Section [43| will need to be augmented 
with topological calculations on the phase boundary, e.g. the Euler characteristic, as 
was carried out in the recent work of Tcramoto and Nishiura [3 6) . 
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